home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / dbesk.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  7.8 KB  |  221 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((nulim (make-array 2 :element-type 'f2cl-lib:integer4)))
  12.   (declare (type (simple-array f2cl-lib:integer4 (2)) nulim))
  13.   (f2cl-lib:fset (f2cl-lib:fref nulim (1) ((1 2))) 35)
  14.   (f2cl-lib:fset (f2cl-lib:fref nulim (2) ((1 2))) 70)
  15.   (defun dbesk (x fnu kode n y nz)
  16.     (declare (type (simple-array double-float (*)) y)
  17.              (type f2cl-lib:integer4 nz n kode)
  18.              (type double-float fnu x))
  19.     (prog ((w (make-array 2 :element-type 'double-float)) (cn 0.0) (dnu 0.0)
  20.            (elim 0.0) (etx 0.0) (flgik 0.0) (fn 0.0) (fnn 0.0) (gln 0.0)
  21.            (gnu 0.0) (rtz 0.0) (s 0.0) (s1 0.0) (s2 0.0) (t_ 0.0) (tm 0.0)
  22.            (trx 0.0) (xlim 0.0) (zn 0.0) (i 0) (j 0) (k 0) (mz 0) (nb 0) (nd 0)
  23.            (nn 0) (nud 0))
  24.       (declare (type f2cl-lib:integer4 nud nn nd nb mz k j i)
  25.                (type (simple-array double-float (2)) w)
  26.                (type double-float zn xlim trx tm t_ s2 s1 s rtz gnu gln fnn fn
  27.                 flgik etx elim dnu cn))
  28.       (setf nn (f2cl-lib:int-sub (f2cl-lib:i1mach 15)))
  29.       (setf elim (* 2.303 (- (* nn (f2cl-lib:d1mach 5)) 3.0)))
  30.       (setf xlim (* (f2cl-lib:d1mach 1) 1000.0))
  31.       (if (or (< kode 1) (> kode 2)) (go label280))
  32.       (if (< fnu 0.0) (go label290))
  33.       (if (<= x 0.0) (go label300))
  34.       (if (< x xlim) (go label320))
  35.       (if (< n 1) (go label310))
  36.       (setf etx
  37.               (coerce (the f2cl-lib:integer4 (f2cl-lib:int-sub kode 1))
  38.                       'double-float))
  39.       (setf nd n)
  40.       (setf nz 0)
  41.       (setf nud (f2cl-lib:int fnu))
  42.       (setf dnu (- fnu nud))
  43.       (setf gnu fnu)
  44.       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
  45.       (setf fn (- (+ fnu n) 1))
  46.       (setf fnn fn)
  47.       (if (< fn 2.0) (go label150))
  48.       (setf zn (/ x fn))
  49.       (if (= zn 0.0) (go label320))
  50.       (setf rtz (f2cl-lib:fsqrt (+ 1.0 (* zn zn))))
  51.       (setf gln (f2cl-lib:flog (/ (+ 1.0 rtz) zn)))
  52.       (setf t_ (+ (* rtz (- 1.0 etx)) (/ etx (+ zn rtz))))
  53.       (setf cn (* (- fn) (- t_ gln)))
  54.       (if (> cn elim) (go label320))
  55.       (if (< nud (f2cl-lib:fref nulim (nn) ((1 2)))) (go label30))
  56.       (if (= nn 1) (go label20))
  57.      label10
  58.       (setf fn gnu)
  59.       (setf zn (/ x fn))
  60.       (setf rtz (f2cl-lib:fsqrt (+ 1.0 (* zn zn))))
  61.       (setf gln (f2cl-lib:flog (/ (+ 1.0 rtz) zn)))
  62.       (setf t_ (+ (* rtz (- 1.0 etx)) (/ etx (+ zn rtz))))
  63.       (setf cn (* (- fn) (- t_ gln)))
  64.      label20
  65.       (if (< cn (- elim)) (go label230))
  66.       (setf flgik -1.0)
  67.       (multiple-value-bind
  68.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7)
  69.           (dasyik x gnu kode flgik rtz cn nn y)
  70.         (declare (ignore var-0 var-1 var-2 var-3 var-6 var-7))
  71.         (setf rtz var-4)
  72.         (setf cn var-5))
  73.       (if (= nn 1) (go label240))
  74.       (setf trx (/ 2.0 x))
  75.       (setf tm (/ (+ gnu gnu 2.0) x))
  76.       (go label130)
  77.      label30
  78.       (if (= kode 2) (go label40))
  79.       (if (> x elim) (go label230))
  80.      label40
  81.       (if (/= dnu 0.0) (go label80))
  82.       (if (= kode 2) (go label50))
  83.       (setf s1 (dbesk0 x))
  84.       (go label60)
  85.      label50
  86.       (setf s1 (dbsk0e x))
  87.      label60
  88.       (if (and (= nud 0) (= nd 1)) (go label120))
  89.       (if (= kode 2) (go label70))
  90.       (setf s2 (dbesk1 x))
  91.       (go label90)
  92.      label70
  93.       (setf s2 (dbsk1e x))
  94.       (go label90)
  95.      label80
  96.       (setf nb 2)
  97.       (if (and (= nud 0) (= nd 1)) (setf nb 1))
  98.       (multiple-value-bind
  99.           (var-0 var-1 var-2 var-3 var-4 var-5)
  100.           (dbsknu x dnu kode nb w nz)
  101.         (declare (ignore var-0 var-1 var-2 var-3 var-4))
  102.         (setf nz var-5))
  103.       (setf s1 (f2cl-lib:fref w (1) ((1 2))))
  104.       (if (= nb 1) (go label120))
  105.       (setf s2 (f2cl-lib:fref w (2) ((1 2))))
  106.      label90
  107.       (setf trx (/ 2.0 x))
  108.       (setf tm (/ (+ dnu dnu 2.0) x))
  109.       (if (= nd 1) (setf nud (f2cl-lib:int-sub nud 1)))
  110.       (if (> nud 0) (go label100))
  111.       (if (> nd 1) (go label120))
  112.       (setf s1 s2)
  113.       (go label120)
  114.      label100
  115.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  116.                     ((> i nud) nil)
  117.         (tagbody
  118.           (setf s s2)
  119.           (setf s2 (+ (* tm s2) s1))
  120.           (setf s1 s)
  121.           (setf tm (+ tm trx))
  122.          label110))
  123.       (if (= nd 1) (setf s1 s2))
  124.      label120
  125.       (f2cl-lib:fset (f2cl-lib:fref y (1) ((1 *))) s1)
  126.       (if (= nd 1) (go label240))
  127.       (f2cl-lib:fset (f2cl-lib:fref y (2) ((1 *))) s2)
  128.      label130
  129.       (if (= nd 2) (go label240))
  130.       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
  131.                     ((> i nd) nil)
  132.         (tagbody
  133.           (f2cl-lib:fset (f2cl-lib:fref y (i) ((1 *)))
  134.                          (+
  135.                           (* tm
  136.                              (f2cl-lib:fref y
  137.                                             ((f2cl-lib:int-sub i 1))
  138.                                             ((1 *))))
  139.                           (f2cl-lib:fref y ((f2cl-lib:int-sub i 2)) ((1 *)))))
  140.           (setf tm (+ tm trx))
  141.          label140))
  142.       (go label240)
  143.      label150
  144.       (if (= kode 2) (go label160))
  145.       (if (> x elim) (go label230))
  146.      label160
  147.       (if (<= fn 1.0) (go label170))
  148.       (if (> (* (- fn) (- (f2cl-lib:flog x) 0.693)) elim) (go label320))
  149.      label170
  150.       (if (= dnu 0.0) (go label180))
  151.       (multiple-value-bind
  152.           (var-0 var-1 var-2 var-3 var-4 var-5)
  153.           (dbsknu x fnu kode nd y mz)
  154.         (declare (ignore var-0 var-1 var-2 var-3 var-4))
  155.         (setf mz var-5))
  156.       (go label240)
  157.      label180
  158.       (setf j nud)
  159.       (if (= j 1) (go label210))
  160.       (setf j (f2cl-lib:int-add j 1))
  161.       (if (= kode 2) (go label190))
  162.       (f2cl-lib:fset (f2cl-lib:fref y (j) ((1 *))) (dbesk0 x))
  163.       (go label200)
  164.      label190
  165.       (f2cl-lib:fset (f2cl-lib:fref y (j) ((1 *))) (dbsk0e x))
  166.      label200
  167.       (if (= nd 1) (go label240))
  168.       (setf j (f2cl-lib:int-add j 1))
  169.      label210
  170.       (if (= kode 2) (go label220))
  171.       (f2cl-lib:fset (f2cl-lib:fref y (j) ((1 *))) (dbesk1 x))
  172.       (go label240)
  173.      label220
  174.       (f2cl-lib:fset (f2cl-lib:fref y (j) ((1 *))) (dbsk1e x))
  175.       (go label240)
  176.      label230
  177.       (setf nud (f2cl-lib:int-add nud 1))
  178.       (setf nd (f2cl-lib:int-sub nd 1))
  179.       (if (= nd 0) (go label240))
  180.       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nd)))
  181.       (setf gnu (+ gnu 1.0))
  182.       (if (< fnn 2.0) (go label230))
  183.       (if (< nud (f2cl-lib:fref nulim (nn) ((1 2)))) (go label230))
  184.       (go label10)
  185.      label240
  186.       (setf nz (f2cl-lib:int-sub n nd))
  187.       (if (= nz 0) (go end_label))
  188.       (if (= nd 0) (go label260))
  189.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  190.                     ((> i nd) nil)
  191.         (tagbody
  192.           (setf j (f2cl-lib:int-add (f2cl-lib:int-sub n i) 1))
  193.           (setf k (f2cl-lib:int-add (f2cl-lib:int-sub nd i) 1))
  194.           (f2cl-lib:fset (f2cl-lib:fref y (j) ((1 *)))
  195.                          (f2cl-lib:fref y (k) ((1 *))))
  196.          label250))
  197.      label260
  198.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  199.                     ((> i nz) nil)
  200.         (tagbody (f2cl-lib:fset (f2cl-lib:fref y (i) ((1 *))) 0.0) label270))
  201.       (go end_label)
  202.      label280
  203.       (xermsg "SLATEC" "DBESK" "SCALING OPTION, KODE, NOT 1 OR 2" 2 1)
  204.       (go end_label)
  205.      label290
  206.       (xermsg "SLATEC" "DBESK" "ORDER, FNU, LESS THAN ZERO" 2 1)
  207.       (go end_label)
  208.      label300
  209.       (xermsg "SLATEC" "DBESK" "X LESS THAN OR EQUAL TO ZERO" 2 1)
  210.       (go end_label)
  211.      label310
  212.       (xermsg "SLATEC" "DBESK" "N LESS THAN ONE" 2 1)
  213.       (go end_label)
  214.      label320
  215.       (xermsg "SLATEC" "DBESK" "OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL" 6
  216.        1)
  217.       (go end_label)
  218.      end_label
  219.       (return (values nil nil nil nil nil nz)))))
  220.  
  221.